/* Get data on birthrates */

import delimited using ca_birth_rate.csv, clear
drop if county == "California"
rename county countystr
encode county, gen(county)

reshape long births, j(year) i(county)

tempfile birthrates
save `birthrates'



joinby county year using `birthrates', unm(using)

drop _merge

joinby county year using fpr_mapdata



/* Table showing pollution results of the STAR experiments*/

eststo all_ca: estpost summarize fpr_orig fpr_ideal_sim  co_pred co_sim_ideal  nox_pred nox_sim_ideal  pm10_pred pm10_sim_ideal  if year == 2009
eststo la_county: estpost summarize fpr_orig fpr_ideal_sim  co_pred co_sim_ideal  nox_pred nox_sim_ideal  pm10_pred pm10_sim_ideal  if year == 2009 & county == 19
eststo bay_area: estpost summarize fpr_orig fpr_ideal_sim  co_pred co_sim_ideal  nox_pred nox_sim_ideal  pm10_pred pm10_sim_ideal  if year == 2009 & inlist(county,1,7,21,28,38,41,43,48,49)
eststo san_diego: estpost summarize fpr_orig fpr_ideal_sim  co_pred co_sim_ideal  nox_pred nox_sim_ideal  pm10_pred pm10_sim_ideal  if year == 2009 & county == 37

#delimit ;
esttab all_ca la_county bay_area san_diego using star_experiment.tex, booktabs replace main(mean) not
varlabels(fpr_orig "\quad Baseline" fpr_ideal_sim `"\quad Re-assign using Quality Score"' 
nox_pred "\quad Baseline" nox_sim_ideal `"\quad Re-assign using Quality Score"' 
pm10_pred "\quad Baseline" pm10_sim_ideal `"\quad Re-assign using Quality Score"' 
co_pred "\quad Baseline" co_sim_ideal `"\quad Re-assign using Quality Score"' , 
blist(fpr_orig `"Quality Score Score\\"' nox_pred "\addlinespace \nox\ Level (PPB)\\" co_pred "\addlinespace CO Level (PPB) \\" pm10_pred "\addlinespace PM\$_{10}\$ Level (\$\mu/m^3\$) \\"))
nonum noobs nostar nonote mtitles("All CA" "LA County" "Bay Area" "San Diego")
title(Predicted 2009 Quality Score and Pollution, Simulating an Ideal and Real STAR Program\label{tab:star-experiment})
note(Notes: Baseline scenario shows actual average Quality Score and predicted CO and \nox\ from the model in column 2 of table \ref{tab:fpr-interaction}.  \`\`Reassign using Quality Score'' redistributes inspections at stations with Quality Score below 0.4 to stations which pass these measures, recalculates the county average Quality Score accordingly, and predicts CO and \nox.  \`\`Reassign to STAR Stations'' redistributes inspections at stations failing the SVFR and FPR threshholds to stations which pass these measures.)
substitute( "\begin{tabular}" "\setlength{\linewidth}{.1cm}\newcommand{\contents}{\begin{tabular}"
            "\end{tabular}" "\end{tabular}}\setbox0=\hbox{\contents}\setlength{\linewidth}{\wd0-2\tabcolsep-.25em}\contents"
            "{l}{\footnotesize" "{p{\linewidth}}{\footnotesize")
;
#delimit cr

/* Calculate costs and benefits*/



/* Costs of Smog Check: $20 per inspection, plus 20 minutes testing time times $25/hr times total number of inspections, plus repair costs*/

gen smogcheck_cost = N*(20 + .33*23.82) + repair_cost*num_repairs

/* Benefits of Smog Check: 17.6 avoided infant deaths per 100,000 births per 1000 units CO (Currie et al 2009), plus 10 avoided infant deaths per 100,000 births per 
unit PM10 (Knittel et al 2015); all times $7.6M VSL as used by EPA*/

gen smogcheck_co_benefits = (17.6 * (births/100000) * (co_sim_nosmog - co_pred)/1000)*7600000
gen smogcheck_pm_benefits = (10 * (births/100000) * (pm10_sim_nosmog - pm10_pred))*7600000
gen smogcheck_benefits = smogcheck_co_benefits + smogcheck_pm_benefits


sort county year
format star* smogcheck* %15.0fc

list county star_* smogcheck* if year == 2009

tabstat star* smogcheck* if year == 2009 & smogcheck_benefits!=., stat(sum) format(%15.0fc) col(stats) varwidth(32)
tabstat star* smogcheck* if year == 2002 & smogcheck_benefits!=., stat(sum) format(%15.0fc) col(stats) varwidth(32)

foreach var in  smogcheck_co_benefits smogcheck_pm_benefits smogcheck_cost{
	replace `var' = `var'/1000000
}

#delimit ;
graph bar (sum) smogcheck_co_benefits smogcheck_pm_benefits smogcheck_cost 
if year >=2002 & year <=2009 & smogcheck_benefits!=.,
over(year) legend(label(1 "CO Benefits") label(2 "PM Benefits") label(3 Costs))
graphregion(color(white))
b1title(Year) ytitle("Social Costs and Benefits (millions)")
;
#delimit cr
graph export no_smogcheck_sim.pdf, replace